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ABSTRACT 

How  to  reduce  the  horizontal  pressure  gradient  error  is  a  key  issue  of  using  (7-coordinate  ocean  models, 
especially  of  using  primitive  equation  models  for  coastal  regions.  The  error  is  caused  by  the  splitting  of  the 
horizontal  pressure  gradient  term  into  two  parts  and  the  subsequent  incomplete  cancellation  of  the  truncation 
errors  of  those  parts.  Due  to  the  fact  that  the  higher  the  order  of  the  difference  scheme,  the  less  the  truncation 
error  and  the  more  complicated  the  computation,  a  sixth-order  difference  scheme  for  the  (7-coordinate  ocean 
models  is  proposed  in  order  to  reduce  error  without  increasing  complexity  of  the  computation.  After  the  analytical 
error  estimation,  the  Semi-spectral  Primitive  Equation  Model  is  used  to  demonstrate  the  benefit  of  using  this 
scheme.  The  stability  and  accuracy  are  compared  with  those  of  the  second-order  and  fourth-order  schemes  in 
a  series  of  calculations  of  unforced  flow  in  the  vicinity  of  an  isolated  seamount.  The  sixth-order  scheme  is 
shown  to  have  error  reductions  by  factors  of  5  compared  to  the  fourth-order  difference  scheme  and  by  factors 
of  50  compared  to  the  second-order  difference  scheme  over  a  wide  range  of  parameter  space  as  well  as  a  great 
parametric  domain  of  numerical  stability. 


1.  Introduction 

In  shallow-water  prediction  models  the  effects  of  bot¬ 
tom  topography  must  be  taken  into  account.  This  can 
be  done  by  using  a  terrain-following  cr-coordinate  sys¬ 
tem,  where  the  water  column  is  divided  into  the  same 
number  of  grid  cells  independent  of  depth.  Let  (x%,  y*, 
z )  denote  Cartesian  coordinates  and  (x,  y,  a)  sigma  co¬ 
ordinates.  In  most  sigma  coordinate  ocean  models  the 
relationship  between  the  two  coordinate  systems  is 


where  g  is  the  acceleration  due  to  gravity  and  p  is  the 
density  of  the  fluid.  Substitution  of  (3)  into  (2)  and 
use  of 

d  2d 


dz  H  da 

yield 

dp  dp  1  dH 

—  =  —  -  -0  -  o)pg—, 
OX*  ox  1  ox 


(4) 


z 

x  =  **,  y  =  y*,  a  =  \  +  2— — -  (1) 

H(x,  y) 

where  z  and  a  increase  vertically  upward  such  that  z  = 
0,  cr  =  1  at  the  surface  and  z  =  ~H  and  a  =  —  1  at 
the  bottom;  H  =  H(x,  y)  is  the  bottom  topography. 

We  restrict  our  attention  to  two  dimensions  (x,  a)  for 
simplicity.  Horizontal  gradients  in  the  z-  and  cr-coor- 
dinates  are  related  by 
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Most  ocean  models  are  hydrostatic  balanced;  that  is. 
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which  is  used  for  computing  the  horizontal  pressure 
gradient. 

A  problem  has  long  been  recognized  in  computing 
the  horizontal  pressure  gradient  in  the  cr-coordinate  sys¬ 
tem  (e.g.,  Gary  1973;  Haney  1991;  Mellor  et  al.  1994; 
McCalpin  1994)  the  horizontal  pressure  gradient  be¬ 
comes  a  difference  between  two  terms,  which  leads  to 
a  large  truncation  error  at  steep  topography.  We  refer 
the  reader  to  a  comprehensive  review  by  Haney  (1991). 
Mellor  et  al.  (1994)  obtained  an  error  estimation  scheme 
and  showed  that  the  pressure  gradient  error  is  advec- 
tively  eliminated  after  a  long  time  integration.  McCalpin 
(1994)  employed  a  fourth-order  difference  scheme  to 
eliminate  the  pressure  gradient  error  by  a  factor  of 
10-20.  Following  McCalpin’s  path,  we  propose  a  sixth- 
order  difference  scheme  for  a  further  reduction  of  pres¬ 
sure  gradient  error. 
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2.  Second-,  fourth-,  and  sixth-order  schemes 

In  ocean  models,  the  most  common  spatial  difference 
scheme  used  to  approximate  the  first  derivative  is  a  two 
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point,  centered,  second-order  staggered  scheme  such  as 
the  C-grid  scheme  (Hadivogel  et  al.  1991)  and  the 
B-grid  scheme  (Robinson  1993): 

/V\  P,+m  ~  Pi— 1/2  _  1  /d3/A  A2 
\dx/.  A  24\dx3). 


Pi  = 


Pi- 1/2  Pi+ 1/2 

2 


+  0(  A2), 


(5) 


where  A  is  the  grid  spacing.  The  truncation  error  of  this 
scheme  is  0( A2). 

For  the  fourth-order  C-grid  scheme. 


Pi- 3/2  ~  27 p,_1/2  +  27pi+m  -  Pi+3/2 

[dxj.  24A 


9  1 
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10  lo 
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the  truncation  error  is  0(A4).  Comparing  (6)  with  (5), 
the  second-order  difference,  the  error  ratio  between  the 
fourth-order  and  second-order  schemes  is  estimated  by 
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Since  the  truncation  error  decreases  with  the  increase 
of  the  order  of  the  difference  scheme,  it  might  be  ben¬ 
eficial  to  use  an  even  higher-order  difference  scheme. 
Thus,  we  propose  a  sixth-order  difference  scheme 


(<>P\ 
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to  compute  the  horizontal  pressure  gradient.  For  the  two 
points  nearest  to  the  boundary,  we  use  the  lower-order 
schemes  (taking  one-side  boundary  as  an  example), 


dp\  ^  Pm  -  Pm 

dx):  A 
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Comparing  (8)  with  (6),  the  error  ratio  between  the 
sixth-order  and  fourth-order  schemes  is  estimated  by 


5  (d7P 
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3.  Seamount  test  case 

a.  Model  description 

Suppose  a  seamount  located  inside  a  periodic /-plane 
(/„  =  1 0  4  s  ')  channel  with  two  solid,  free-slip  bound¬ 
aries  along  constant  y.  Unforced  flow  over  a  seamount 
in  the  presence  of  resting,  level  isopycnals  is  an  ideal 
test  case  for  the  assessment  of  pressure  gradient  errors 
in  simulating  stratified  flow  over  topography.  The  flow 
is  assumed  to  be  reentrant  (periodic)  in  the  alongchannel 
coordinate  (i.e.,  x  axis).  We  use  this  seamount  case  of 
the  Semi-spectral  Primitive  Equation  Model  (SPEM) 
version  3.9  to  test  the  new  difference  scheme.  The  reader 
is  referred  to  the  original  reference  (Haidvogel  et  al. 
1991)  and  the  SPEM  3.9  User’s  Manuel  (Fledstrom 
1995)  for  datailed  information.  In  the  horizontal  direc¬ 
tions  the  model  uses  the  C-grid  and  the  second-order 
finite  difference  discretization  except  for  the  horizontal 
pressure  gradient  for  which  the  user  has  a  choice  of 
either  second-order  or  fourth-order  difference  discreti¬ 
zation  (McCalpin  1994).  In  the  vertical  direction  the 
model  uses  a  boundary-fitted  cr-coordinate  system.  The 
discretization  is  by  spectral  collocation  using  Cheby- 
shev  polynomials.  Our  model  configuration  is  similar 
to  that  of  Beckmann  and  Flaidvogel  (1993)  and  Mc¬ 
Calpin  (1994).  The  time  step  and  grid  size  used  here 
are 
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Fig.  1.  The  seamount  geometry. 


At  =  675  s. 


Ax  =  Ay  =  5  km. 


b.  Topography 

The  domain  is  a  periodic  channel,  320  km  long  and 
320  km  wide.  The  channel  walls  are  solid  (no  normal 
flow)  with  free-slip  viscous  boundary  conditions.  The 
channel  has  a  far-field  depth  /zmax  and  in  the  center  in¬ 
cludes  an  isolated  Gaussian-shape  seamount  with  a 
width  L  and  an  amplitude  hs  (Fig.  1), 


h(x,  y)  =  h„ 


hsex  p 


(x  -x0)2  +  (y  -  y0): 


,  (10) 


where  (x0,  y0)  are  the  longitude  and  latitude  of  the  sea¬ 
mount  center.  The  far-field  depth  (/zmax)  is  fixed  as  5000 
m.  But  the  seamount  amplitude  (hs)  changes  from  500 
to  4500  m,  and  the  lateral  scale  of  the  seamount  (L) 
varies  from  6  to  40  km  for  the  study. 


c.  Density  field 

Suppose  that  the  background  fluid  is  at  rest  and  with 
an  exponentially  stratified  mean  density 

p(z)  =  28-2  exp >  (H) 

where  z  is  the  vertical  coordinate,  and  Hp  =  1000  m. 
Here  a  constant  density,  1000  kg  m~3,  has  been  sub¬ 
tracted  for  error  reduction.  When  a  density  disturbance 
p\  =  p  exp (z/Hp)  was  introduced  initially, 

p,  =  m  +  pi,  (i2) 


the  fluid  was  slightly  less  stably  stratified  than  the  ref¬ 
erence  field  for  each  computation.  Here  p  is  called  the 
density  anomaly.  The  larger  the  p ,  the  less  stable  the 
fluid  will  be.  In  this  study,  p  varies  from  0.1  to  1  kg  m  3. 

Following  Beckmann  and  Haidvogel  (1993)  and 
McCalpin  (1994),  we  subtract  the  mean  density  field 
p(z)  before  integrating  the  density  field  to  obtain  pres¬ 
sure  from  the  hydrostatic  equation. 


d.  Lateral  viscosity 

Ideally,  the  new  difference  scheme  should  be  tested 
with  no  lateral  diffusion  of  density.  This  is  due  to  the 
fact  that  the  density  diffusion  along  cr  surfaces  generates 
horizontal  gradients  wherever  the  cr  surfaces  are  not  flat, 
and  then  produces  horizontal  pressure  gradients  that 
drive  currents  in  much  the  same  way  as  the  pressure 
gradient  errors  (McCalpin  1994).  Unfortunately,  the  ab¬ 
sence  of  the  horizontal  diffusion  keeps  the  small-scale 
pressure  disturbances  generated  by  topographic-scale 
density  advection  and  the  induced  small-scale  velocity 
fields,  which  in  turn  cause  the  computational  instability 
problem.  Thus,  some  lateral  viscosity  on  cr  surfaces  is 
required  in  the  momentum  equations  to  maintain  sta¬ 
bility.  A  constant  coefficient  ( AH )  biharmonic  formu¬ 
lation  is  used  here  for  the  lateral  viscosity,  which  varies 
from  1010  to  10s  m4  s~‘  in  this  study. 

4.  Standard  case 

a.  Model  parameters 

In  this  section,  we  assess  the  sixth-order  scheme  (8) 
by  the  cumulative  effects  of  the  pressure  gradient  errors. 
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Fig.  2.  Instantaneous  error  pattern  after  5  days  of  integration  for  the  unforced  (exponential  stratification)  experiment  with  sixth-order 
difference  scheme:  (a)  u,  (b)  v,  (c)  w,  and  (d)  mass  transport  streamfunction.  Here  u,  v,  w  are  evaluated  at  the  slice  (facing  upchannel) 
through  the  center  of  the  seamount. 


A  variety  of  configuration  were  used  for  the  test.  These 
will  be  described  in  next  section.  At  first,  we  set  up  a 
standard  test  case  to  compare  errors  among  three  dif¬ 
ferent  (second,  fourth,  and  sixth)  schemes: 

L  =  40  km,  hs  =  4500  m, 

p  =  0.2  kg  m~3,  Ah  =  1010  m4  s-1.  (13) 

The  SPEM3.9  standard  was  integrated  for  20  days  for 


the  standard  test  case  using  ordinary  and  compact 
fourth-  and  sixth-order  difference  schemes  for  comput¬ 
ing  horizontal  pressure  gradient.  Figure  2  displays  errors 
in  the  streamfunction  and  velocity  fields  after  perform¬ 
ing  5  days  of  integration  using  the  sixth-order  scheme. 
The  mass  transport  streamfunction  has  a  large-scale 
eight-lobe  pattern  centered  on  the  seamount.  This  sym¬ 
metric  structure  can  be  found  in  all  the  fields.  After  5 
days  of  integration,  the  model  generates  spurious  cur¬ 
rents  of  0(0.001  cm  s->). 
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Fig.  3.  Peak  error  velocity  for  the  second-,  fourth-,  and  sixth-order 
schemes. 


b.  Temporal  variations  of  peak  error  velocity 

Owing  to  a  very  large  number  of  calculations  per¬ 
formed,  we  discuss  the  results  exclusively  in  terms  of 
the  maximum  absolute  value  of  the  spurious  velocity 
(called  peak  error  velocity)  generated  by  the  pressure 
gradient  errors.  Figure  3  shows  the  time  evolution  of 
the  peak  error  velocity  for  the  first  20  days  of  integration 
with  the  second-,  fourth-,  and  sixth-order  schemes.  The 
peak  error  velocity  fluctuates  rapidly  during  the  first  few 
days  integration.  After  the  5  days  of  integration,  the 
peak  error  velocity  shows  the  decaying  inertial  oscil¬ 
lation  superimposed  on  mean  values:  0.4  cm  s  1  (the 
second  order),  0.042  cm  s  (the  fourth  order),  and  0.01 
cm  s_1  (the  sixth  order).  The  steady  approach  of  the 
peak  error  velocities  to  these  values  for  the  three 
schemes  indicates  the  stability  of  the  computation.  Fur¬ 
thermore,  the  time-mean  peak  error  velocity  for  the  sec¬ 
ond-order  case  is  roughly  40-50  times  that  of  the  sixth- 
order  case. 

c.  Performance  of  the  second-,  fourth-,  and 
sixth-order  schemes 

Feasibility  of  using  the  sixth-order  scheme  is  twofold: 
1)  drastic  error  reduction  and  2)  no  drastic  CPU  time 
increase.  We  use  four  quantities  to  compare  the  errors: 
the  maximum  peak  error  velocity  during  the  10-day  in¬ 
tegration  (V,),  the  averaged  peak  error  velocity  during 
the  5-10  day  integration  ( V2 ),  the  maximum  peak  error 
pressure  gradient  during  the  10-day  integration  (PG,), 
and  averaged  peak  error  pressure  gradient  during  the  5- 
10  day  integration  (PG2).  Here  V,  (or  PG,)  indicates  the 
maximum  error  during  the  geostrophic  adjustment  stage, 
and  V2  (or  PG2)  shows  the  mean  error  after  the  adjust¬ 
ment.  Usually,  V,  >  V2,  and  PG,  >  PG2.  The  perfor- 


Table  1.  Ratios  of  CPU  time  and  various  errors  between  the 
second-order  and  the  fourth  (sixth)-order  schemes. 


Ratio 

CPU 

u 

U 

PG, 

PG, 

2nd/4th 

0.9020 

8.873 

10.111 

11.682 

14.306 

2nd/6th 

0.9087 

44.643 

46.083 

76.336 

81.301 

mance  of  each  scheme  is  listed  in  Table  1.  Both  V ,  and 
V2  reduce  their  magnitude  by  factors  around  10  from 
the  second-order  scheme  to  the  fourth-order  scheme  (V, 
from  7.4  X  10_1  to  8.17  X  10~2  cm  s_1,  V,  from  4.61 
X  10_1  to  4.43  X  10~2  cm  s_1X  and  by  factors  around 
5  from  the  fourth-order  scheme  to  the  sixth-order 
scheme  (V,  from  8.17  X  10~2  to  1.61  X  10~2  cm  s_1, 
V2  from  4.43  X  10~2  to  9.75  X  10~3  cm  s-1)-  The 
pressure  gradient  error  reduces  even  more  when  the 
high-order  schemes  are  used:  PG,  reduces  from  9.7  X 
10  1  to  8.3  X  10  2  N  m  3  from  the  second-order 
scheme  to  the  fourth-order  scheme  (a  factor  of  11.7) 
and  from  8.3  X  10~2  to  1.27  X  10~2  N  m~3  from  the 
fourth-order  scheme  to  the  sixth-order  scheme  (a  factor 
of  6.5).  PG2  reduces  from  7.1  X  10~3  to  4.96  X  10~4 
N  m  3  from  the  second-order  scheme  to  the  fourth- 
order  scheme  (a  factor  of  14.3)  and  from  4.96  X  1CU4 
to  8.75  X  1CU5  N  m~3  from  the  fourth-order  scheme 
to  the  sixth-order  scheme  (a  factor  of  5.7).  The  CPU 
time  increases  10%  from  the  second-order  to  fourth- 
order  and  sixth-order  scheme.  Notice  that  there  is  no 
CPU  increase  from  the  fourth-order  to  the  sixth-order 
scheme,  which  shows  the  great  potential  of  using  the 
sixth-order  scheme. 

5.  Sensitivity  studies 

We  performed  various  sensitivity  study  to  show  the 
dependence  of  model  results  on  the  four  parameters  list¬ 
ed  in  (13).  For  each  sensitivity  study,  we  vary  only  one 
parameter  and  then  integrate  the  model  for  10  days  using 
the  second-,  fourth-,  and  sixth-order  schemes,  respec¬ 
tively. 

a.  Seamount  width 

Seamount  widths  (L)  of  40,  20,  15,  12,  and  10  km 
were  used.  A  summary  of  V ,  and  V2  after  10  days  of 
integration  for  each  case  is  presented  in  Fig.  4.  The 
seamount  with  L  larger  than  10  km  can  be  resolved  by 
the  grid  spacing  used  here  (5  km).  The  figure  clearly 
shows  a  great  improvement  of  using  high-order 
schemes.  This  improvement  enhances  as  L  increases.  At 
L  =  10  km,  V,  reduces  by  a  factor  of  1.45  from  the 
second-order  scheme  to  the  fourth-order  scheme  and  a 
factor  of  1.18  from  the  fourth-order  scheme  to  the  sixth- 
order  scheme.  At  L  =  40  km,  V,  reduces  by  a  factor  of 
9.06  from  the  second-order  scheme  to  the  fourth-order 
scheme  and  a  factor  of  5.08  from  the  fourth-order 
scheme  to  the  sixth-order  scheme.  Similar  reduction  of 
V2  was  found  (Fig.  4)  from  the  second-order  scheme  to 
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Fig.  4.  Sensitivity  of  error  on  the  seamount  horizontal  scale  for  the 
second-,  fourth-,  and  sixth-order  schemes:  (a)  Vx  and  (b)  V2. 


the  fourth-order  scheme,  and  from  the  fourth-order 
scheme  to  the  sixth-order  scheme. 

b.  Seamount  amplitude 

Seamount  amplitudes  (hs)  of  4500,  4000,  3500,  3000, 
2000,  1000,  and  500  m  were  used.  A  summary  of  Vx 
and  V2  after  10  days  of  integration  for  each  case  is 
presented  in  Fig.  5.  The  figure  clearly  shows  a  great 
improvement  using  high-order  schemes.  This  improve¬ 
ment  enhances  as  hs  increases.  At  hs  =  500  m  (i.e.,  a 
small  seamount  with  a  slope  of  0.0125),  reduces  by 
a  factor  of  2.55  from  the  second-order  scheme  to  the 
fourth-order  scheme  and  a  factor  of  1.10  from  the 
fourth-order  scheme  to  the  sixth-order  scheme.  At  hs  = 
4500  m  (i.e.,  a  large  seamount  with  a  slope  of  0.1125), 
V,  reduces  by  a  factor  of  9.06  from  the  second-order 
scheme  to  the  fourth-order  scheme  and  of  5.08  from  the 
fourth-order  scheme  to  the  sixth-order  scheme.  Similar 
reduction  of  V2  was  found  (Fig.  5)  from  the  second- 


Fig.  5.  Sensitivity  of  error  on  the  seamount  amplitude  for  the  sec¬ 
ond-,  fourth-,  and  sixth-order  schemes:  (a)  V,  and  (b)  V2. 


order  scheme  to  the  fourth-order  scheme,  and  from  the 
fourth-order  scheme  to  the  sixth-order  scheme. 

c.  Density  anomaly 

Density  anomalies  (p)  of  1,  0.8,  0.6,  0.4,  0.2,  and  0.1 
kg  m  3  were  used.  A  summary  of  V,  and  V2  after  10 
days  of  integration  for  each  case  is  presented  in  Fig.  6. 
The  figure  clearly  shows  a  great  improvement  using 
high-order  schemes.  This  improvement  is  not  sensitive 
to  the  change  of  p.  At  p  =  0.1  kg  m~3  (i.e.,  a  small 
density  anomaly),  reduces  by  a  factor  of  9.05  from 
the  second-order  scheme  to  the  fourth-order  scheme  and 
a  factor  of  5.08  from  the  fourth-order  scheme  to  the 
sixth-order  scheme.  At  p  =  1  kg  m~3  (i.e.,  a  large  density 
anomaly),  V ,  reduces  by  a  factor  of  9.10  from  the  sec¬ 
ond-order  scheme  to  the  fourth-order  scheme  and  of 
5.07  from  the  fourth-order  scheme  to  the  sixth-order 
scheme.  Similar  reduction  of  V2  was  found  (Fig.  6)  from 
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Fig.  6.  Sensitivity  of  error  on  the  density  anomaly  for  the  second-, 
fourth-,  and  sixth-order  schemes:  (a)  V,  and  (b)  V2- 


the  second-order  scheme  to  the  fourth-order  scheme,  and 
from  the  fourth-order  scheme  to  the  sixth-order  scheme. 

d.  Horizontal  viscosity 

Horizontal  viscosities  ( AH )  of  1010,  1 09,  and  10s  m4 
s  1  were  used.  A  summary  of  V,  and  V2  after  10  days 
of  integration  for  each  case  is  presented  in  Tables  2  and 
3,  which  clearly  shows  a  great  improvement  using  high- 
order  schemes.  This  improvement  is  not  very  sensitive 
to  the  change  of  AH.  At  AH  =  1010  m4  s-1,  V,  reduces 
by  a  factor  of  9.06  from  the  second-order  scheme  to  the 
fourth-order  scheme  and  a  factor  of  5.08  from  the 
fourth-order  scheme  to  the  sixth-order  scheme.  At  AH 
=  108  m4  s  *,  V,  reduces  by  a  factor  of  9.19  from  the 
second-order  scheme  to  the  fourth-order  scheme  and  of 
5.10  from  the  fourth-order  scheme  to  the  sixth-order 
scheme.  Similar  reduction  of  V2  was  found  (Table  3) 
from  the  second-order  scheme  to  the  fourth-order 
scheme,  and  from  the  fourth-order  scheme  to  the  sixth- 
order  scheme. 


Table  2.  Sensitivity  of  V{  (cm  s  *)  on  AH  for  the  second-,  fourth-, 
and  sixth-order  schemes. 


Ah  (m4  s-) 

Scheme 

10s 

IO5 

1010 

2nd 

8.994  X  101 

8.696  X  10  1 

7.399  X 

10-' 

4th 

9.782  X  10“2 

9.527  X  10-2 

8.171  X 

io-2 

6th 

1.919  X  HU2 

1.869  X  10-2 

1.610  X 

io-2 

Table  3. 

Sensitivity  of  V2  (cm  s  *)  on  AH  for  the  second-, 
and  sixth-order  schemes. 

fourth-, 

Ah  (m4  s-) 

Scheme 

10s 

IO5 

1010 

2nd 

5.364  X  10-' 

5.120  X  10  ' 

4.613  X 

10-' 

4th 

4.819  X  10“2 

4.681  X  10-2 

4.431  X 

10-2 

6th 

9.333  x  10“3 

8.503  X  10-3 

9.748  X 

10-3 

6.  Conclusions 

1)  The  cr-coordinate,  pressure  gradient  error  depends 
on  the  choice  of  difference  schemes.  By  choosing  an 
optimal  scheme,  we  may  reduce  the  error  a  great  deal 
without  increasing  the  horizontal  resolution.  Analytical 
analysis  shows  that  the  truncation  error  of  the  fourth- 
order  scheme  may  be  1-2  orders  of  magnitude  smaller 
than  the  second-order  scheme,  and  the  truncation  error 
of  the  sixth-order  scheme  may  be  1-2  orders  of  mag¬ 
nitude  smaller  than  the  fourth-order  scheme. 

2)  The  SPEM  Version  3.9  is  used  to  demonstrate  the 
benefit  of  using  the  sixth-order  scheme.  A  series  of  cal¬ 
culations  of  unforced  flow  in  the  vicinity  of  an  isolated 
seamount  are  performed.  The  results  show  that  the  sixth- 
order  scheme  has  error  reductions  by  factors  of  5  com¬ 
pared  to  the  fourth-order  difference  scheme,  and  by  fac¬ 
tors  of  50  compared  to  the  second-order  difference 
scheme  over  a  wide  range  of  parameter  space  as  well 
as  a  great  parametric  domain  of  numerical  stability. 

3)  Using  the  sixth-order  scheme  does  not  require 
much  more  CPU  time.  Taking  SPEM3.9  as  an  example, 
the  CPU  time  for  the  sixth-order  scheme  is  almost  the 
same  as  for  the  fourth-order  scheme,  and  10%  more 
than  for  the  second-order  scheme. 

4)  Since  the  fourth-order  different  scheme  has  error 
reductions  by  factors  of  10  compared  to  the  second- 
order  difference  scheme,  there  is  no  real  advantage  to 
going  to  a  higher-order  scheme  if  the  bottom  topography 
is  not  too  complicated.  The  need  for  more  accuracy  will 
go  up  with  increasingly  complex  bottom  topography  on 
small  scales,  so  one  might  expect  that  future  demand 
for  accuracy  will  increase  as  models  strive  for  more 
realism. 
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